#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <omp.h>

void met_gauss_jacob(double **,double *,int ,int ,int ,double);

int main(int argc,char *argv[]){

        int J_ORDER, J_ITE_MAX, J_ROW_TEST;  /**Variaveis do programa**/
        double J_ERROR;

        int i,j;

        double **MA;
        double *VB;

        scanf("%d",&J_ORDER);    /**Ordem da matriz**/
        scanf("%d",&J_ROW_TEST); /**Local da linha a ser realizado o teste da linha**/
        scanf("%lf",&J_ERROR);   /**erro permitido**/
        scanf("%d",&J_ITE_MAX);  /**Numero máximo de interações**/

        MA = (double **)malloc(sizeof(double *)*J_ORDER); /**Aloca a matriz A[][]**/
        for(i = 0;i < J_ORDER;i++)
                MA[i] = (double *)malloc(sizeof(double)*J_ORDER);

        VB = (double *)malloc(sizeof(double)*J_ORDER);    /**Aloca o vetor B[]**/

         #pragma omp for
        for(i = 0;i < J_ORDER;i++){  /**Coloca os dados do arquivo na matriz A**/
                for(j = 0;j < J_ORDER;j++)
                        scanf("%lf",&MA[i][j]);
        }

         #pragma omp for
        for(i = 0;i < J_ORDER;i++)/**Coloca os dados do arquivo no vetor B**/
                scanf("%lf",&VB[i]);


        met_gauss_jacob(MA,VB,J_ORDER,J_ITE_MAX,J_ROW_TEST,J_ERROR);


        free(MA);
        free(VB);

        return(0);
}

int criterio_linha(double **M,int J_ORDER){
        int i,j;
        double sum = 0.0,max = 0.0;
        //  omp_set_num_threads(2);

        #pragma omp for reduction(+=sum) /**Verifica a condição do método**/
        for(i = 0;i < J_ORDER;i++){
                sum = 0.0;
                for(j = 0;j < J_ORDER;j++){
                        if(i != j)
                                sum += M[i][j];
                }
                sum = sum/M[i][i];
                if(max < sum)
                        max = sum;
        }

        if(max < 1)
                return 1;
        else
                return 0;
}

void met_gauss_jacob(double **A,double *B,int J_ORDER,int J_ITE_MAX,int J_ROW_TEST,double J_ERROR){
        double *X;
        double *X1;
        //double *Dif;
         // omp_set_num_threads(2);
        double numerator = 0.0, denominator = 0.0;
        int count = 0,i,j;

        double div,dif,result=0.0;  

        for(i = 0;i < J_ORDER;i++)
                X = (double *)malloc(sizeof(double)*J_ORDER);

        for(i = 0;i < J_ORDER;i++)
                X1 = (double *)malloc(sizeof(double)*J_ORDER);

        for(i = 0;i < J_ORDER;i++) 
                X[i] = 0;

        if(criterio_linha(A,J_ORDER)){ /**Verifica o critério do método**/

                do{

                        count++;
                        denominator = 0.0;
                        numerator = 0.0;
                        #pragma omp for
                        for(i = 0;i < J_ORDER;i++){ /**Realiza a multiplicação da Matriz A pelo vetor B**/
                                        X1[i] = 0.0;
                                for(j = 0;j < J_ORDER;j++){
                                        if(i != j)
                                                X1[i] = X1[i] + A[i][j] * X[j];
                                }
                                X1[i] = 1 / A[i][i] * (B[i] - X1[i]);

                                dif = X1[i] - X[i];

                                if(dif < 0){
                                        dif *= (-1);
                                }

                                if(numerator < dif)
                                        numerator = dif;

                                dif = X1[i];

                                if(dif < 0)
                                        dif *= (-1);

                                if(denominator < dif)
                                        denominator = dif;
                        }
                        #pragma omp for
                        for(i = 0;i < J_ORDER;i++)/**Atribui o X(k+1) em X(k)**/
                                X[i] = X1[i];

                        div = numerator/denominator;

                }while(div >= J_ERROR && count < J_ITE_MAX);

            #pragma omp for
                for(i = 0; i < J_ORDER; i++){
             result += A[J_ROW_TEST][i]*X1[i]; /**Acumula os resultados**/
        }
        printf("Iterations: %d\n", count); /**Printa os resultados**/
        printf("RowTest: %d => [%lf] =? %f\n", J_ROW_TEST, result, B[J_ROW_TEST]);
        }
        else{
                printf("Não obedece aos criterios.\n");
        }

        free(X);
        free(X1);
}
